Compaction and dilation rate dependence of stresses in gas-fluidized beds 
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A particle dynamics-based hybrid model, consisting of monodisperse spherical solid particles and 
volume-averaged gas hydrodynamics, is used to study traveling planar waves (one-dimensional trav- 
eling waves) of voids formed in gas-fluidized beds of narrow cross sectional areas. Through ensemble- 
averaging in a co-traveling frame, we compute solid phase continuum variables (local volume fraction, 
average velocity, stress tensor, and granular temperature) across the waves, and examine the rela- 
tions among them. We probe the consistency between such computationally obtained relations and 
constitutive models in the kinetic theory for granular materials which are widely used in the two- 
fluid modeling approach to fluidized beds. We demonstrate that solid phase continuum variables 
exhibit appreciable "path dependence" , which is not captured by the commonly used kinetic theory- 
based models. We show that this path dependence is associated with the large rates of dilation and 
compaction that occur in the wave. We also examine the relations among solid phase continuum 
variables in beds of cohesive particles, which yield the same path dependence. Our results both for 
beds of cohesive and non-cohesive particles suggest that path-dependent constitutive models need 
to be developed. 

PACS numbers: 



I. INTRODUCTION 

Flows involving gas-particle mixtures are ubiquitous in 
nature and in engineering practice 0, 0, 0, 01 • Homoge- 
neous flows are often unstable, and spatiotemporal struc- 
tures such as traveling waves, clusters and streamers of 
particles, and bubble-like voids are commonly observed 
in fluidized beds. The origin of the instabilities leading 
to these structures has been studied extensively in the 
literature [100011. 

The number of particles present in these flows is large, 
rendering detailed description of the motion of all the 
particles impractical. Hence, macroscopic flow charac- 
teristics are often probed through analysis of continuum 
models derived by averaging the equations governing 
the motion of the individual particles and the intersti- 
tial fluid. In this approach, the fluid and solid phases 
are treated as inter-penetrating continua, and locally- 
averaged quantities such as the volume fractions and ve- 
locities of fluid and solid phases appear as dependent 
variables. The details of flow at the level of individual 
particles appear in the averaged equations through the 
effective fluid and particle phase stresses and the inter- 
phase interaction force, for which one must postulate con- 
stitutive models. It is now known that inhomogeneous 
structures such as traveling waves and bubble-like voids 
in fluidized beds can be reproduced in a qualitatively 
correct manner through very simple constitutive models 
for the above-mentioned quantities (see e.g., Glasser et 
al. ^). Quantitative predictions are still elusive in many 
cases. 
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The derivation of constitutive models for these terms 
has been studied extensively in the literature. For exam- 
ple, the kinetic theory of granular materials (KTGM) for 
the quantiflcation of the stresses in a rapidly shearing as- 
sembly of particles undergoing only binary collisions has 
been discussed by many authors [1 H HlgllTollTlL IT^ . 
and it has found widespread applications in many gas- 
particle flow problems 0, . Many of the early stud- 
ies of the kinetic theory focused on spherical, uniformly 
sized, non-cohesive particles, while more recently there 
have been attempts to extend it to sli ghtl y cohesive sys- 
tems (see e.g., Kim and Arastoopour [lj|)- Particle dy- 
namics simulations of shear flow of nearly homogeneous 
assemblies of non-cohesive particles have been used to 
validate the kinetic theory [ll El El [13 ■ 

As mentioned above, gas-particle flows in fluidized 
beds readily form inhomogeneous structures, where the 
particle assembly undergoes compaction and dilation in 
an alternating periodic manner while also undergoing a 
shear. For example, in fluidized beds, this results from 
the flow patterns created by a succession of bubble-like 
voids rising through the bed. Quite frequently, such ef- 
fects occur over length scales in the range of 10 - 50 par- 
ticle diameters [1 E3 , and the effect of the rate of com- 
paction or dilation accompanying inhomogeneous flows 
at such length scales on the rheological characteristics is 
largely unknown. 

While the KTGM does include a bulk viscosity term 
to account for the effect of compaction/dilation rate on 
the solid phase pressure, the adequacy of this correction 
remains untested. A good understanding of the typical 
rates of compaction and dilation that can occur in flu- 
idized gas-particle suspensions and their influence on the 
stresses is a necessary step in the quest for constitutive 
models that can lead to quantitative predictions of their 
flow behavior. The flrst goal of this study is to ana- 
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lyze a simple, idealized gas-particle flow problem where 
a particle assembly undergoes compaction and dilation 
in an alternating manner and also manifests spatial in- 
homogeneity that is typical of fluidized suspensions, and 
probe how the stresses are influenced by compaction or 
dilation events. To this end, we have carried out particle 
dynamics-based hybrid simulations of a periodic assem- 
bly of uniformly sized, spherical, non-cohesive, frictional 
and inelastic particles fluidized by a gas, and created an 
inhomogeneous flow pattern that takes the form of a one- 
dimensional (vertically) traveling wave (ID-TW). In this 
approach, solid phase stresses are determined through 
particle-particle interactions, rather than by postulated 
constitutive models, and various inter-particle forces such 
as cohesive forces can readily be accounted for with- 
out introducing additional assumptions. By traveling 
along this wave and averaging over many realizations, we 
have extracted the spatial variation of various ensemble- 
averaged quantities. It will be shown that the stresses 
manifest strong path dependence. When the assembly un- 
dergoes dilation, the stress levels are very small and are 
essentially independent of the rate of dilation; in contrast, 
when compaction occurs, the stresses depend strongly on 
the rate of compaction. 

We have also extended these simulations to include 
inter-particle attractive forces; specifically, we considered 
a simple model to capture van der Waals force between 
the particles. It will be shown through an analysis of 
traveling waves in modestly cohesive systems that the 
path dependence is even greater in cohesive systems. 

The rest of the paper is organized as follows. The 
model we use, a particle dynamics model coupled with 
volume-averaged gas hydrodynamics, is described in 
Sec.m The traveling wave structures and properties ob- 
tained in our simulations of beds of non-cohesive particles 
are presented in Sec. lIII Al In Sec. lIIIBl solid phase con- 
tinuum variables across the waves are analyzed and as- 
sumptions used in KTGM are probed. The waves formed 
in beds of moderately cohesive particles are analyzed in 
Sec. nil CI followed by conclusions in Sec. IIVI 



II. METHOD 

A. DEM-based hybrid model 

The solid phase is modeled as a collection of discrete 
particles of "soft" monodisperse spheres (the deforma- 
tion is accounted for by overlaps), whose individual tra- 
jectories are computed by integrating classical equations 
of motion. Such an approach is often referred to as the 
soft-sphere molecular dynamics simulation or the discrete 
element method (DEM). When a pair of objects (par- 
ticles or system boundaries) get into contact with each 
other, the interaction force is analyzed in the normal and 
tangential directions separately, Fcont = (F„,Fs). We 
adopt the so-called spring-dashpot model of Cundall and 



Strack using a Hookean spring 

F„ = (fcnA„ - 7„'y„)n, (1) 
Fs = -sign{vs) X mm{ktAs,fJ.\Fn\)s, (2) 

where kn is the spring stiffness in the normal direction 
(n is the corresponding unit vector, pointing from the 
contact point toward the particle center); A„ is the over- 
lap of the particles in the normal direction; 7„ is the 
damping coefficient determined by fc„ and the normal 
coefficient of restitution e characterizing the inelastic ki- 
netic energy loss upon contact; and v„ is the normal com- 
ponent of the relative velocity at contact. The interac- 
tion in the tangential direction is modeled by a spring 
and a slider, where the tangential stiffness kt is deter- 
mined by kn and the Poisson's ratio of the material i^p 
[a value of 0.3 is used for all the computations in our 
study; kt — 2A;„(1 — vp)/{2 — vp)]; Ag is the tangential 
displacement from the initial contact; = Vg • s is the 
tangential component of the relative velocity at contact 
[vs = n X (vy X n)]; and s is the unit vector in the tan- 
gent plane coUinear with the component of the relative 
velocity at contact. The total tangential force is limited 
by the Coulomb frictional force /i|r„|, where /i is the co- 
efficient of friction. More details on the model equations 
can be found elsewhere p^ . 

The equation of motion for each particle accounts for 
the gas phase effects following the volume-averaged hy- 
drodynamics miilllil: 

'^p-^ = + Fcont + Fc + -^/5(0) (Ug - Vp) - VpVp, 

(3) 

where nip and Vp are individual particle mass and ve- 
locity, respectively. The right-hand side includes various 
forces acting on the particle, the last two terms arising 
from the gas-solid two-way coupling; the total force act- 
ing on the particles due to the ffuid is commonly parti- 
tioned into the local drag part and the effective buoyant 
part, as was done here (see e.g.. Ye et al. [l^). The 
first term is the body force due to gravity, where g is the 
gravitational acceleration, and the second term is afore- 
mentioned contact force. The third term accounts for 
the van der Waals cohesive force Fc ~ — ^^fi, foUow- 
ing Hamaker's formula |2^ in the limit oi s <^ ri — r2 
(monodisperse spheres), where A is the Hamaker con- 
stant, s is the inter-surface distance, and ri and r2 are the 
(identical) radii of particles undergoing interaction. We 
avoid the singularity of this formula at contact by intro- 
ducing a minimum separation distance 5* of 0.4 nm J^, 
below which the same attractive force is applied; i.e., 
Fc = max[— ^p-, — ^-p^p-Jn. The fourth term accounts 
for the drag force, where (3 is the interphase momentum 
transfer coefficient, (f) is the local particle phase volume 
fraction, is the local-average gas phase velocity, and Vp 
is the individual particle volume. The last term accounts 
for the hydrodynamic force due to the gradually varying 
part of the pressure field, where p is the local-average gas 
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phase pressure. 

In general, the gas phase variables in the last two terms 
are obtained by simultaneously integrating the balance 
equations for gas-solid mixtures; however, for the case of 
ID incompressible gas phase, Eq. can be significantly 
reduced by considering the ID continuity relation: 

(l-0)ug + 0u, = U,, (4) 

and the reduced momentum balance equation: 

O = -(l-0)Vp + /3(0)(u,-Ug), (5) 

where Uj is the coarse-grained particle velocity, and Ug 
is the superficial gas flow velocity. After some manipula- 
tion [using Eqs. I@J and jSJ], Eq. Q can be reduced to 
the following one, where gas phase effects appear as addi- 
tional terms in the individual equation of motion, involv- 
ing solid phase coarse-grained, continuum variables p^ : 



dt 



V, 



(6) 



For the momentum transfer coefficient /3, we use an 
empirical expression proposed by Wen and Yu p^ : 



4 dp 



2.65 



(7) 



where C^i is the drag coefficient, Pg is the gas phase mass 
density, and dp is the particle diameter. The drag coeffi- 
cient proposed by Rowe [23| is employed 



Cd 
where 



24 ( 

0.44, 



0.f5i?e0-687) 



Reg < 1000, 
Reg > 1000, 



(1 - (l))Pgdp\Ug - Us 



(8) 



(9) 



and Pg is the gas phase viscosity. We further simplify 
these formulas using the assumption of Reg <SC I, as we 
consider only small particles. 



For the fully fluidized states that we are interested in, 
the total weight of the bed is supported by the pressure 
drop, which can be written as 



P\z=0 - P\z=L = Psg4>avgL 



(10) 



where ps is the mass density of the solid phase, g = |g| is 
the acceleration due to gravity, and 4'avg is the average 
volume fraction of the bed. Solving Eq. IjlUI) together 
with the ID continuity relation Eq. and the momen- 
tum balance equation Eq. jsj yields an appropriate value 
for the superficial gas flow velocity Us that fully fluidizes 
the bed in this small periodic box. We use this value of 
Us in Eq. ®. 



Nondimensionalization 



Casting Eqs. © and in a dimensionless form, us- 
ing Ps, dp, \fgdp, \J dpi g as characteristic density, length, 
velocity, and time, one obtains the following nondimen- 
sional groups (arrows indicate changes in the notation 
from dimensional variables to nondimensional variables 
that will be used henceforth): 







spring stiffness 


Psgdl ' 

Us 


Us <- 


superficial gas flow rate 


y/gdp 


(5 = 


5* 
dp' 

A 


scaled minimum separation distance 


Bo = 
St = 


cohesive Bond number 
Stokes number 


A-Kpsgd^S^ ' 
Psg Up 

Pg 



together with nondimensional parameters, namely e, p, 
Up, and L the nondimensional bed height. We character- 
ize particles by usual values (e = 0.9; p = 0.1; vp = 0.3) 
in experiments, except when we consider "ideal" parti- 
cles in Figs. 151 and If 81 In the following sections, all the 
results are presented in nondimensional form. 



D. Computation of solid phase continuum variables 



B. Fully periodic box 

Traveling waves in a bed of a finite depth (as in exper- 
imental systems) are hardly perfectly periodic. In order 
to eliminate such an imperfection from the wave, we con- 
sider an idealized geometry of a fully periodic (in all three 
directions) small box of height L which is commensurate 
with the wavelength A in a deep fiuidized bed The 
periodic boundary conditions for both lateral directions 
are achieved with the usual geometrical wrapping. In the 
direction of gravity, we adopt the following procedure: 



Solid phase coarse-grained continuum variables at ID 
discrete grid points separated by are computed by 
distributing the particle mass and momenta values to two 
nearby grid points using a halo function h that linearly 
decreases to zero around the ith particle located at z = 



h{z] Zj) 



I - |z-z,|/Az 




for \z - Zi\ < Az, 
otherwise, ^ ' 



where Az is the grid spacing. It is readily seen that h 
has the property that the solid phase quantities of each 
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FIG. 1: A schematic of a halo function in Eq. IIH . ap- 
plied to a particle (a large circle) whose center is located at 
z = Zi, which is used to compute solid phase coarse-grained 
continuum variables. Three dots along the abscissa represent 
nearby grid points. The halo function linearly decreases from 
the unity aX z = Zi to zero a.t z — Zi ± Az. 



particle are distributed to the two nearby grid points, 
inversely proportional to the distance to each grid point. 
The number density n[= {6/Tr)(f>/dp] and u^, on the grid 
point zq are then defined simply as 



N 



n{za) = ^h{zo;zi), (12) 

1=1 

N 

i{zo)us{za) = ^h{zo;zi)vp^i, (13) 




0.5 0.5 



FIG. 2: (Color online) Left panel: A snapshot of non-cohesive 
particles {Bo = 0) in a fully periodic (in all three directions) 
box (10 X 10 X 18) whose height is commensurate with one 
wavelength of a one-dimensional traveling wave (left panel). 
Right panel: Corresponding coarse-grained solid phase con- 
tinuum variables: local volume fraction (p, vertical component 
of the locally averaged velocity , granular temperature T, 
and the average normal stress Tr(^)/_D, where D = 3 is the 

dimension. [All quantities are dimensionless.] 



Az/10 throughout the bed. 



where Zi is the ith particle location. Following the same 
procedure, it is straightforward to compute the granular 
temperature tensor T: 



Z = ((Vp,i - Us) (8) (Vp,i - Us)) , 



(14) 



or the scalar granular temperature T [— j^Tt{X)], where 

D is the dimensionality, (E" is the dyadic tensor product, 
and "Tr" denotes the trace (of a tensor) . 

The full stress tensor consists of a kinetic or dynamic 
part and a virial or static part (2^ : 



1 

V 



^ ® 1, 



cev 



(15) 



where \p^i = Vp^i — Ug is the fluctuating velocity, fc is 
the force between contacting particles 1 and 2, and Ic — 
ri — r2 is the displacement vector between the centers 
of the particles under consideration. The second term is 
summed over all the contacts in the averaging volume V. 

Solid phase continuum variables are computed, us- 
ing the above halo function, on the nodes separated by 
the grid spacing Az. In order to compute smoothly 
varying continuum variables across a bed, this proce- 
dure is repeated on many subgrid points in between 
the grid points, by translating all the nodes by a small 
amount. As a result, all the continuum variables are 
computed on uniformly distributed points separated by 



III. SIMULATION RESULTS 

We consider gas-fluidized beds of non-cohesive as well 
as cohesive particles in fully periodic boxes of a narrow 
square-shaped 10 x 10 cross sectional area. The height 
of the box L, which is the same as the wavelength A in 
a deep fluidized bed, is set to be 18 in most cases; it 
is varied when we study the effect of the wavelength in 
Fig.CI We also set St at 55 in most cases, except when St 
is varied in Fig. |31 Grid spacing for continuum variables 
Az is set to be 1.5 in all the computations we will present. 
We have checked that the quantitative results slightly 
vary when a different value is used for Az, but the main 
results remain the same. We also have confirmed that the 
results are virtually independent of the cross sectional 
area [l9|| . 

We obtain smooth wave profiles by transforming the 
bed to the co-traveling frame (with the wave), and av- 
eraging continuum variables over hundreds of snapshots 
of fully-developed waves at different instances (as well 
as spatial smoothing through the halo function). When 
beds of cohesive particles are considered fSec. IIII C|) . we 
characterize the level of cohesion by the cohesive Bond 
number Bo, which is defined to be the ratio of the maxi- 
mum cohesive force (at contact) to the gravitational force 
acting on the particle, namely the weight of the particle. 
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FIG. 3: (Color online) Continuum variable profiles across 
fully developed waves formed in beds of "realistic" particles 
(solid lines: e — 0.9, = 0.1) and of "ideal" particles (dashed 
lines: e — 1.0, /i — 0.0), where (pavg is set to be 0.44 for 
both cases; there are only slight differences between the two 
cases. The wave travels from left to right, and gravity is acting 
toward the negative z-direction. {Bo — 0; St = 55; fc„ — 
2.0 X 10^.) [All quantities are dimensionless.] 



FIG. 4: (Color online) Continuum variable profiles across 
fully developed waves formed in beds of two types of particles 
with different stiffness (solid lines: fc„ = 2.0 x 10^; dashed 
lines: Kn ~ 2.0 X 10*), where (j)avg is set to be 0.44 for both 
cases; only slight differences between the two cases. The in- 
tegration time stepsize for the case of k„ — 2.0 x 10^ is larger 
by a factor of 30. (e = 0.9; /i = 0.1; Bo = 0; St = 55; <pavg = 
0.44.) [All quantities are dimensionless.] 



A. Traveling planar waves 

When beds are subject to a superficial gas flow veloc- 
ity above the minimum fluidization rate (Us > Um/), 
homogeneously expanded states are often unstable, and 
beds form inhomogeneous structures, such as bubbles 
In beds of narrow cross sectional areas, bubbling occurs 
in the form of propagating planar waves (or ID-TW) of 
voidage p9| , which are reproduced in our model. A snap- 
shot of a bed of non-cohesive particles, forming ID-TW, 
illustrates that particles in the upper plug "rain down" 
through the void region, and accumulate at the top of 
the lower plug. As a result, the void region propagates 
up (Fig. 121 red particles are moving up and blue ones are 
moving down, in online color figure). We use the afore- 
mentioned halo function (spatial smoothing) and time 
averaging, to compute the corresponding continuum vari- 
ables, which are shown in the right panel. For notational 
convenience, we drop a subscript of to represent the 
coarse-grained solid phase velocity by u (and represent 
its z-component by Uz). 

We start by varying the particle dissipation parame- 
ters (e and /i) to probe the sensitivity of the wave pro- 
files to them. Our simulations reveal that the wave 
profiles hardly change with the dissipation parameters 
(Fig. 13). As the particles become more frictional, the 
plateau volume fraction values in the plug get smaller, 
since frictional particles pack more loosely in the lower 
plug because of frictional resistance to rotations (see an 
inset in Fig. However, the change is very small, be- 
cause the inertia of the particles is much larger than 
that of the fluid in gas-fluidized beds (hence particles 
accumulate with larger velocities, compared to those 
in liquid-fluidized beds). It is expected that the effect 
would be more pronounced for smaller particles. Li and 



Kuipers [30|,|3l| showed that heterogeneous structures do 
exist in "ideal" particles (e = 1.0, /i — 0.0), and that dis- 
sipative collisions dramatically intensify the heterogene- 
ity of the structures. We do not see such dramatically in- 
creased effects; the waveform changes only slightly, which 
suggests that the wave mainly arises from the instability 
associated with the interplay between solid phase inertia 
and the nonlinear volume-fraction-dependent gas drag, 
as argued by Jackson (S^ . Inelastic collisions themselves 
indeed lead to inhomogeneous structures through a well- 
known "clustering instability" fsSj; however, our simula- 
tion suggests that this is a secondary effect in the forma- 
tion of ID-TW. 
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FIG. 5: (Color online) Continuum variable profiles across 
fully developed waves (of the same wavelength), formed 
in beds of non-cohesive particles {Bo — 0;e — 0.9; fi = 
0.1;(/f>aua = 0.36) at different Stokes numbers {St = 55 for 
solid lines; 28 for dashed lines; 14 for dot-dashed lines). Vol- 
ume fraction profiles are nearly the same for the three cases, 
while all the other variables decrease in magnitude, as St de- 
creases. [All quantities are dimensionless.] 
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FIG. 6: (Color online) Profiles of solid phase continuum variables for fully developed waves formed in beds of non-cohesive 
particles {Bo = 0;St — 55; e = 0.9; ^ = 0.1) in a fully periodic box of 10 x 10 x 18, shown for three cases of different average 
volume fractions (a) <j>avg = 0.36; (b) (j>avg = 0.44; (c) <j>avg = 0.52. Vertical bars in (a) represent the standard deviation of time- 
averaged data (over 500 snapshots). Different diagonal components of the stress tensor axx (dot-dashed line) and a^z (dashed 
line) (and those of the granular temperature tensor in inset) are shown separately in (b). [All quantities are dimensionless.] 



One ad hoc value used in our simulation is the spring 
stiffness kn (as is the case in most other DEM simula- 
tions). In principle, it can be computed from the Young's 
modulus of the material under consideration; however, 
such values limit the integration step to extremely small 
sizes, which leads to expensive computations. We used 
two values of kn differing by three orders of magnitude 
(while all the other parameters are kept the same), and 
compared the solid phase continuum variables (Fig. 0J. 
We find that the results are nearly the same, except that 
the average normal stress Tr (£)/£> is somewhat larger 

when a stiffer spring constant value is used; however, the 
difference in the stress for the two cases is smaller than 
its level of fluctuation [which is shown later in Fig.|Sl(a)]. 
For the parameters used in Fig. 0] the duration of head- 
on collision Ate = tt/ ^ (^)(1 + -^), where m* is the re- 
duced mas and a — tt/ loge, corresponds to 8.2 x lO^^sec 
and 2.6 x 10"^ for fc„ = 2.0 x 10^ and 2.0 x 10*, respec- 
tively; it is reasonably small even for the smaller value 
of kn- For all the rest of the computations, we will use 
the smaller fc„ value (~ 10^), as it allows us to use an 
integration timestep size 10~^Atc, which scales with 
kn^^"^ as shown above) larger by a factor of 30. 

For the wave profiles shown in Figs. |3| through |7| 1121 
and 1141 the gravity is acting leftward (c/. right panel 
in Fig. 121). In sustained traveling waves, particles in the 
upper plug rain down through the void region, and the 
waves travel to right. In the rest of the paper, we will re- 
fer the region where particles fall off from the bottom of 
upper plug (characterized by duz/dz > 0) as the dilation 
region, and where the particles accumulate at the top of 
lower plug {duz/dz < 0) as the compaction region. Let us 
track the particle phase continuum variables along the di- 
rection of gravity, starting from the dilation region: Both 
the granular temperature and the particle phase stresses 
are negligibly small across the void region. As a stream 
of particles collide with the lower plug (in the compaction 
region), both T and the stress rapidly increase. The par- 
ticles quickly become solid-like in the compaction region 



through coUisional dissipation, and T starts to decrease 
(before the stress does). The collisional stress becomes 
dominant, and it keeps increasing until the volume frac- 
tion reaches its plateau value. Inside a plug the stress 
monotonically decreases, as particles gradually lose con- 
tacts and eventually fall off. 

Before analyzing the wave profiles any further, we com- 
pare them at different Stokes numbers by varying the gas 
phase viscosity (Fig. [S]). For a fixed wavelength, the vol- 
ume fraction profiles change only slightly; as St decreases 
(by increasing by factors of 2 and 4), the plateau value 
decreases and the overall profile gets smoothed a little 
bit. More pronounced gas drag makes the magnitude of 
all other solid phase quantities decrease by comparable 
factors. Quantitative changes in variables (volume frac- 
tion, granular temperature, and average normal stress) 
are better shown in insets of Fig. [S| We will set St = 55 
for the rest of the computations, which corresponds to 
simulating with usual air. 

The volume fraction profiles for a range of average 
volume fraction (j)avg [Fig. El (a) through (c)] illustrate 
that the plateau values in the close-packed regions re- 
main nearly the same, while both the depth (wave am- 
plitude) and the width of the void region gradually de- 
crease with increasing (l)avg (when the wavelength is kept 
the same) . Both the solid phase average normal stress in 
the compaction region (at the top of the lower plug) and 
the maximum average speed, decrease as (pavg increases. 
It is worth noting that the solid phase pressure exhibits 
large fluctuations, whereas all the other variables (0, u, 
and T) exhibit relatively small fluctuations [Fig. (a)]. 
In all cases, in the compaction region, the zz-component 
of the stress tensor is larger than the other components 
{xx- and yy-) which are the same due to symmetry [it 
is shown only in Fig. El (b)]. The observed fluctuations 
of the stress are quantitatively comparable with the dif- 
ference between azz and a^x, however, the difference is 
observed consistently in average quantities for all cases. 
The anisotropy is significantly large in the granular tem- 
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FIG. 7: (Color online) Continuum variable profiles across 
fully developed waves of different wavelength A (A = 12 for 
dot-dashed lines; 18 for solid lines; 24 for dashed lines), 
formed in beds of non-cohesive particles (e = 0.9; fi — 
0.1; Bo — 0; St = 55; (jtavg ~ 0.44). The plateau volume frac- 
tions remain the same, while the amplitudes of the volume 
fraction profiles decrease as A decreases. Vertical thin solid 
hues at 2 = 12 and 18 represent the system boundaries for 
the cases of the same wavelengths. [All quantities are dimen- 
sionless.] 



perature tensor T [see an inset in Fig. El (b)], but it is 

limited to a narrow region in space. The first normal 
stress difference in sfieared granular flows is shown to be 
a Burnett order effect js^l- However, the normal stress 
difference that we observe in a bubbling fluidized bed has 
a different origin (see Appendix for details). 



The continuum variables for waves of different wave- 
lengths, for a fixed value of average volume fraction, are 
shown in Fig.[7| Once the wavelength is sufhciently long, 
the same plateau volume fraction values are obtained for 
all three cases. Solid lines in Fig.[7|correspond to the data 
in Fig. El (b). From Figs. El and [7| we see that the plateau 
volume fraction values are the same for ranges of average 
volume fractions and wavelengths. This is precisely what 
Glasscr et al. J\ obtained in their simulations of contin- 
uum models. In their analysis, as the wavelength was 
decreased, the amplitude of the wave (i.e. of the volume 
fraction) decreased to zero and the uniform state became 
stable at a Hopf bifurcation point. We simulate the waves 
for a range of wavelengths, and compute amplitudes of 
volume fraction profiles (Fig.lSJ. Our analysis yields simi- 
lar results to those which were obtained in the continuum 
model, but the amplitudes at small wavelengths are finite 
and the bifurcations are imperfect because of the discrete 
nature of the model. No further resolution of the data 
at larger wave numbers is available, as the wavelength 
increases discretely by Az(= 1.5) (i.e. A = 4.5 and 6.0 
were used). 



FIG. 8: The wave amplitudes (of the volume fraction profiles) 
obtained for various wavelengths (4.5 < A < 24) at three 
different average volume fractions, {(fiavg ~ 0.54, 0.56, and 
0.58; e = 0.9; fi = 0.1; Bo ^ 0; St ^ 55.) 



B. Particle phase stresses and viscosities, and path 
dependence 

We use the data in Fig. El to computationally obtain 
constitutive relations, and assess consistency with the as- 
sumptions used in the KTGM, which is widely used for 
two- fluid models 0, . 

We start by plotting the scalar granular temperature 
T and average normal stress -^Tr(£.) against the volume 

fraction for the three cases in Fig. El [see Figs. El (a) and 
(b); in (b), the volume fraction in the plateau region re- 
mains the same for a range of solid phase stresses, and 
the exact distinction between dilation and compaction 
branches in this region is elusive]. The computed values 
in the dilation region are nearly the same for all three 
cases; however, those in the compaction region strongly 
depend on (fiavg- Furthermore, they are significantly dif- 
ferent from (and much larger than) those in the dilation 
region. It is clear that both T and -^Tr(£) lie on two dis- 
tinct branches, depending on where the measurements of 
these variables are made. Similar "path dependence" was 
also observed in lattice-Boltzmann simulations of liquid- 
fluidized beds [s^. 

In the previous section, we have shown that the gran- 
ular temperature tensor manifests significant anisotropy, 
which has to be accounted for, in principle, in our anal- 
ysis. However, an appreciable anisotropy is seen only 
in small portion of the compaction region [see inset of 
Fig. El (b)], while the anisotropy [Fig. El (b)] and path 
dependence [Fig. El (b)] in the stress tensor are manifest 
over a wider region. Thus, it is reasonable to believe 
that neither the anisotropy nor the path dependence of 
the stress arise from the anisotropy of the granular tem- 
perature. Therefore, in what follows, we consider only 
the scalar granular temperature. 

For ID flows: 

^Tt{s) ^ Ps - fJ^b^ , (16) 
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FIG. 9: (Color online) (a) The granular temperature, (b) average normal stress, and (c) average normal stress divided by 
the solid volume fraction and the temperature (i.e. the unity plus inelastic dense gas correction in the equation of state) for 
the three cases in Fig. |Sl(»,o : (jj^vg = 0.36; A, v : 't'avg = 0.44; *, x : cpavg = 0.52). Measurements from the compaction 
region where particles accumulate (characterized by duz/dz < 0) are distinguished from those from the dilation region where 
particles rain down [du^/dz > 0); •,A,*'s are obtained from the compaction regions, and 0,^7, x's are obtained from the 
dilation regions. The same symbols will be consistently used through Fig. 1131 and all quantities are dimensionless. 



where Ps is the dimensionless solid phase pressure, and /if, 
is the dimensionless bulk viscosity. According to KTGM, 



=0r[l + 2(l + e)05o('^)], 



(17) 



where go{(j)) is the isotropic "equilibrium" radial distri- 
bution function (RDF) evaluated at contact. We expect 
that a part of the path dependence of -^Tr(£) arises from 

the path dependence of T. We examine whether the first 
term in Eq. (|16|l alone (i.e. neglecting the bulk viscosity 
effect) is enough to explain the path dependence of the 
stress. We rescale the average normal (dimensionless) 
stress by (j)T in Fig. |51(c), which corresponds to the in- 
elastic and dense gas correction (plus unity) in KTGM, 
assuming that the relation in Eq. H17|) holds. We find that 
the path dependence still persists. Therefore, Eq. H17|) 
cannot adequately explain the path dependence. 

Physically, the path dependence of the average normal 
stress can be rationalized as follows: As a granular as- 
sembly gets dilated or compacted, the RDF at contact 
departs from its equilibrium value, and depends on the 
rate of dilation/compaction. When the relevant dimen- 
sionless group, /VT, is negligibly small, the assembly 
equilibrates quickly and remains close to the equilibrium 
configuration. (Note that, in this study, all the results are 
presented in terms of dimensionless quantities. Since we 
use the particle diameter as the characteristic length, this 
dimensionless group becomes ^j^/VT when Uz, z, and T 
are all dimensionless. In terms of dimensional variables, 
this group would be dp^j^/^/T.) However, when this 
quantity is somewhat larger, the RDF at contact is ex- 
pressed as a perturbation from its equilibrium value, and 
the bulk viscosity term accounts for the departure from 
the equilibrium value. [In such a perturbation approach, 
go in Eq. H17() will be the equilibrium value and is only a 
function of particle volume fraction.] Therefore, there is 
some rationale in trying to explain the path dependence 
as the effect of the bulk viscosity. 



We ascertain that the actual RDF at contact is indeed 
path-dependent by directly computing it from our sim- 
ulations. The isotropic RDF at distance r = r„ is |36l| : 



girn) 



N. 



shell 



^^al,g47^^2A^' 



(18) 



where r„ = (n — 1/2) Ar, Uavg is the average number 
density, and N shell is the number of particles lying in a 
thin spherical shell bounded by (n — 1) Ar and nAr. We 
compute the first five values of the RDF (1 < n < 5, 
using Ar = 10~^), and extrapolate them to estimate its 
value at contact r = 0.5 (Fig. IIUI) : for the sake of better 
comparison with the results in Fig. IHI (c) , a prefactor of 
2(1 -I- e)0 is multiphed and unity is added. The results 
vary a little bit depending on the choice of Ar and the 
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FIG. 10: (Color online) Computationally obtained isotropic 
RDF at contact (see text for details), for the three cases in 
Fig. 1^1 which are appropriately scaled for a better comparison 
with Fig. |U](c). The RDFs indeed exhibit path dependence. 
The dashed line is a scaled plot of an analytical formula used 
by Bagnold, go{(j)) = - (<?f>/<?f>max)^'^^], where cpmax = 0.63 
(~ the plateau volume fraction value) was used. 
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FIG. 11: (Color online) Solid phase dimensionless shear vis- 
cosity scaled by dimensionless \/r, obtained for the three dif- 
ferent values of (fjavg in Fig. |^ Insets include the viscosities 
shown in linear- logarithmic scale. 



FIG. 12: (Color online) A dimensionless group relevant to 
the rate of dilation or compaction, ^jf- /Vt, for the cases in 
Fig. |S| (dashed lines: (jtavg = 0.36; solid lines: (jtavg ~ 0.44; 
dot-dashed lines: (pavg ~ 0.52), plotted together with the 
scaled volume fraction profile. [Here z, Uz, and T are all 
dimensionless.] 



extrapolation method. However, in all cases, the path 
dependence persists. Thus we can indeed expect a non- 
negligible bulk viscosity correction. Before discussing the 
bulk viscosity effect, let us look at the shear viscosity. 

We numerically estimate the first derivative of the 
average velocity fields through the central differencing 
scheme, and compute the effective (dimensionless) shear 
viscosity /is from different components of the stress ten- 
sor: 



1 4 duz 

a..--Tr(^) = --M.^ 



(19) 



As the shear viscosity in KTGM depends on vT, we plot 
/is after rescaling it by ^/T (Fig. [TT|l . We find that scaled 
/is's computed from different collapse nearly onto 

common functional forms, but the path dependence is 
very clear. In the KTGM, the shear viscosity has no 
path dependence 0: 



/is = 



50F 



48(1 
4 



ej.go 



e)<^.go 



r.9o(l 



(20) 



which holds, again, for small mag nitude of^/VT, and 
the apparent path dependence again signals clearly a 
marked departure from conditions assumed in the de- 
velopment of KTGM. 

We have computed the dimensionless quantity ^^/VT 
from our simulation results; see Fig. El Note that the 
magnitude of this quantity is not small compared to 
unity, and so the assumption that the actual RDF at con- 
tact can be written as a small perturbation from the equi- 
librium value, which is routinely made in the theories, 
does not hold in this rather simple inhomogeneous flow. 
The path dependence of the shear viscosity is thus clearly 
associated with effects which are nonlinear in ^j^/Vt. 

Having found that is path-dependent, we can ex- 



pect the bulk viscosity to differ in the compaction 
and dilation branches as well. Strictly speaking, there is 
no unique way of partitioning -^Tr(£.) into Ps and /ife^j^, 

if fib is going to depend on (j>, T, and the path. However, 
it is clear from Fig.^that T and j)Tr{g_} in the dilation 

branch are nearly identical for the three cases, and inde- 
pendent of I \/T. This suggests that the bulk viscosity 
in the dilation branch is close to zero, so that ^Tr(£) m 

the dilation branch can be taken as simply equal to ps 
in that branch. As seen in Fig. E] the shear viscosity in 
the dilation branch is much smaller than that in the com- 
paction branch. Thus, it is not unreasonable to suspect 
that the bulk viscosity correction in the dilation branch 
is small: 

Ps\dil ~ -^Tr(CT)|rfi;, 

where the subscript dil indicates evaluation on the dila- 
tion branch. Armed with this, we can estimate the bulk 
viscosity in the compaction branch in at least two differ- 
ent ways. First, we can set (assume) Ps\dii ~ Ps\comp and 
write 



jj Tr(^) I cojnp 



-^Tr(£)|d»/ 



duz 
dz 



(21) 



where the subscript comp indicates evaluation on the 
compaction branch. Bulk viscosity values in the com- 
paction branch estimated in this manner are shown in 
Fig. El (a) . Alternately one can recognize that the gran- 
ular temperature is different in the two branches; ac- 
counting for the difference, we can write 



1 



1 



-^Tr(CT)|co,„p = -^Tr(£)|dj/ x ' 



T 



T, 



dil 



dUz 

dz 



(22) 

Figure El (b) shows the bulk viscosity values estimated 
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FIG. 14: (Color online) Profiles of solid phase continuum 
variables obtained for non-cohesive {Bo = 0, solid lines) and 
slightly cohesive (Bo = 2, dashed lines) particles, with the 
same average volume fraction, (j)avg = 0.44. The inset shows 
a blow-up of the region where the stress of the cohesive bed 
becomes tensile, (e = 0.9; n = 0.1; St = 55.) 



Effect of cohesion 



FIG. 13: (Color online) Solid phase dimensionless bulk vis- 
cosity scaled by dimensionless VT, obtained for the three dif- 
ferent values of <l)avg in Fig. |31 estimated via (a) Eq. 12H , and 
(b) Eq. 1221 . Insets include the viscosities shown in linear- 
logarithmic scale. 



via Eq. In both Figs. [T^ (a) a nd (b), the bulk vis- 

cosity values are scaled by ^jTcomp- The insets in these 
figures show the results in a linear-logarithmic scale. The 
second approach leads to smaller estimates for the bulk 
viscosity. However, it is clear that irrespective of the 
approach (assumption) taken [i.e. Eq. or if^ ]. 

the estimated values of bulk viscosity in the compaction 
branch are found to be appreciably larger than the shear 
viscosity (compare Figs. II II and IT!^. It is interesting to 
contrast this with KTGM predictions 4]: 

^^b^\4>^9o{l + e)\- (23) 

3 V TT 

which is smaller than the magnitude of the shear viscos- 
ity. 

It should be emphasized that the ID-TW studied here 
is not a peculiar problem representing extreme condi- 
tions. Such sharp volume fraction gradients routinely 
occur around bubble-like voids and clusters in fluidized 
beds, turbulent beds, and fast fluidized beds 0. In such 
devices, the particle assembly is frequently subjected to 
local dilation and compaction in alternating manner. In 
these flows, dpV • u/^/T may not be small and this can 
lead to the type of effects seen in this study. 



It is well known that beds of flne powders are difS- 
cult to fluidize because of the cohesive force between the 
particles. Beds of cohesive powders have attracted in- 
creased attention in recent years; however, their theo- 
retical understanding is still limited. In this section, we 
vary the level of cohesion to probe how the cohesive force 
changes the wave profiles and the relations among con- 
tinuum variables. 

The waves formed in beds of weakly cohesive particles 
{Bo <~ 2) retain a virtually invariant shape over time, 
and travel with a well-defined constant speed. They be- 
come more and more irregular and slow down, until they 
vanish at Bo 8, depending on A and 4'avg', the propa- 
gation speed at Bo = 6 (3) is about 60% (85%) of that in 
beds of Bo = The wave profiles for non-cohesive 
and weakly cohesive particles {Bo < 2) are shown in 
(Fig. I14|l . There are two major differences between the 
waves in beds of non-cohesive and of cohesive particles: 
Firstly, the plateau volume fraction is smaller in cohe- 
sive beds, as cohesive particles pack more loosely than 
non-cohesive ones, due to reasons similar to those for 
frictional particles (See e.g., Dong et al. 1^3 )• Even for 
slightly cohesive particles, this effect is appreciable (see 
inset of Fig. I14() : i.e. the effect is more pronounced com- 
pared to that of a changing friction coefficient (Fig. |3J|. 
Secondly, the average normal stress becomes tensile and 
remains so for a certain range of the height in the dila- 
tion region. At the bottom of the upper plug, the average 
normal stress monotonically decreases until it reaches its 
minimum value (the largest tensile stress) and the assem- 
bly breaks up (inset in Fig. The magnitude of the 
tensile stress cannot exceed the tensile strength of the 
cohesive assembly; thus, one can expect that the former 
corresponds to the tensile strength of this fluidized bed. 
When we estimate the tensile strength via this reason- 
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FIG. 15: (Color online) Solid phase granular temperature and average normal stress for three cases of different cohesion level 
(•, o : Bo = 0; A,v : Bo = 1; *, x : Bo = 2); (pavg = 0.44 in all cases (e = 0.9, fi = 0.1). •, A, *'s are obtained from the 
compaction regions, and o, ^j, x's are obtained from the dilation regions. {St = 55.) The wave speed monotonically slows down 
with increasing level of cohesion, until the wave disappears at Bo ~ 8; the wave speed in beds of Bo — 2 is ^ 90% of that in 
Bo = 0. 



ing, it is almost an order of magnitude smaller that what 
is predicted by Rumpf's model This makes sense, 

since the cohesive assembly in a fluidized bed breaks up 
through the direction of the weakest linkage, in contrast 
to all directions isotropically, as assumed in Rumpf's 
model js^l- We observe further decrease (increase) of 
the plateau volume fraction (the magnitude of maximum 
tensile strength), as Bo increases, until the waves vanish 
(i.e. the bed does not get fluidized any more) at ~ 8. 
The size of agglomerates broken off of the cohesive as- 
sembly however has a distribution, and the wave profile 
becomes irregular for Bo >~ 3, which makes an quanti- 
tative analysis difficult. For sufficiently cohesive particles 
(Bo >~ 8), as one can readily expect, the maximum ten- 
sile stress driven solely by gravity fails to reach the tensile 
strength of the material, and the bed does not form a sta- 
ble traveling wave any more. At that point, an additional 
mechanism, such as mechanical vibration or fine powder 
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coating, has to be implemented to fluidize them |4Ct l41l|. 

As in the beds of non-cohesive particles, both T and 
Tt{s.)/D form lobes when plotted against </>, and manifest 

the same path dependence (Figs . [T^ and [T7|l . Clearly, our 
arguments on pronounced bulk viscosity effects, and on 
path dependence of both the bulk and shear viscosities 
inferred for beds of non-cohesive particles should remain 
valid for beds of moderately cohesive particles as well. 
As Bo increases (in a narrow range of Bo < 2, at a fixed 
average volume fraction), both the maximum granular 
temperature [Fig. 1151 (a)] and the maximum volume frac- 
tion in the plateau region [Fig. El (b)] decrease notice- 
ably, but the maximum (compressive) value attained by 
the average normal stress remains nearly the same. We 
compute the shear viscosity for the three cases shown in 
Fig. El and find that it increases as Bo increases (only 
the values on compaction branches are shown in Fig. 1161 
for clarity), as one can expect. 

For a fixed level of cohesion {Bo — 2), the maximum 
granular temperature decreases as the average volume 
fraction increases [Fig. El (a)] , because the void region 
gets narrower and particles accumulate more gently as in 
the beds of non-cohesive particles with varying average 
volume fraction. The average normal stress value for dif- 



ferent 



-'avg 



in the compaction branch varies. However, 



t I , t • • J 
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the granular temperature [Fig. 1171 (a)] and average nor- 
mal stress [Fig.EI(b)] in the dilation branch for the three 
cases follow (roughly) the same curve until the assembly 
breaks off; as in beds of non-cohesive particles, dilation 
branches follow nearly the same curves. 



FIG. 16; (Color online) Solid phase dimensionless shear vis- 
cosity in the compaction branch {(pavg = 0.44), scaled by di- 
mensionless VT, obtained from the three cases of different 
values of Bo in Fig.[TSl(» : Bo = 0; A : Bo = 1; * ; Bo = 2.) 
The viscosities are shown in linear-logarithmic scale in the 
inset. 



IV. DISCUSSION 

Soft sphere models for particle-particle interactions 
have been used in the literature for nearly three decades 
to analyze a variety of dense particulate flow prob- 
lems [IE 111 il 111 El El 111 111 El El. The 
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FIG. 17: (Color online) Solid phase granular temperature and average normal stress in beds of cohesive particles {Bo = 2; e = 
0.9; /i = 0.1) for three different average volume fractions (•, o : (pavg ~ 0.40; A,v : (pavg ~ 0.44; *, x : 4'avg ~ 0.48). •, A, *'s 
are obtained from the compaction regions, and o, •\/, x's are obtained from the dilation regions. {St = 55.) 



particle dynamics-based hybrid model employed in our 
study, which combines a soft-sphere model for the parti- 
cles with local-average equations for the gas phase, has 
been used previously to simulate flows o f ga s-particle 
mixtures [IS El 121 IH |H3, El El El El, in- 
cluding those occurring in multi-dimensional fluidized 
beds ll m El El E3, El El, and successfully com- 
pared against experimental data El El|- Thus there 
is ample basis for using this approach to probe the na- 
ture of inhomogeneous gas-particle flows. We have used 
this model to simulate one of the simplest spatiotem- 
poral structures possible in gas-fluidized beds, namely 
one-dimensional traveling waves (ID-TW). We believe 
that such a simple structure is ideal for our study, as 
it contains both compacting and dilating regions, and its 
steady nature (in the co-traveling frame) allows us to per- 
form suitable ensemble-averaging, which is essential for 
the analysis of results of discrete simulations. Through 
an analysis of "computational data" generated through 
simulations of these traveling wave structures, we have 
examined in this study the possible consequences of com- 
paction and dilation rates on the stresses generated in the 
particle assembly. Such volume changes accompany most 
gas-particle flows, and understanding their effects on the 
stresses is thus of practical importance. 

It is well established in the literature, through stability 
analysis based on continuum two-fluid models, that such 
traveling waves are the most dominant modes through 
which homogeneously fluidized suspensions lose stability. 
In real gas-fluidized beds, such traveling waves are not 
observed, as they quickly give way to a secondary insta- 
bility leading to bubble-like voids d H E3- However, 
one can readily generate such traveling waves in compu- 
tations by restricting the cross-section of the bed to only 
several particle diameters, thereby suppressing the sec- 
ondary instability H^. Although the ID-TW generated 
in this manner are artificial structures, they represent a 
simple, suitable inhomogeneous flow pattern where the 
particle assembly undergoes compaction and dilation al- 
ternately in a periodic manner. 



We began by generating computational data on ID- 
TW for assemblies of uniformly sized, spherical, non- 
cohesive particles. We demonstrated that the particle 
volume fraction profiles observed in the hybrid model 
calculations are qualitatively similar to those seen pre- 
viously in analyses of two- fluid models ^7^ i5T^. These 
traveling waves in dense fluidized beds take the form of 
vertically rising plugs of particles. The particle volume 
fraction in the plugs was found to be essentially inde- 
pendent of the average particle volume fraction in the 
wave, while the height of the void and the lowest volume 
fraction attained in this void depended systematically on 
it (see Figs. and [T)) — just as in the case of similar 
structures obtained via two-fluid models supplemented 
with ad hoc closures 0, E3 ■ The present hybrid model 
replaces the ad hoc closures for solid phase stresses in the 
two-fluid model with particle-particle interaction models. 

The variation of wave amplitude with wavenumber 
(Fig. |HJ) obtained in our simulations is also similar to that 
obtained with two-fluid models. Typically, in two-fluid 
models, the homogeneous state is stable for short wave- 
lengths, and it gives way to traveling waves via a Hopf bi- 
furcation at some critical wavelength. As the wavelength 
increases further, the wave amplitude increases Al- 
though the discrete simulations undertaken in the present 
study do not exhibit a perfect bifurcation from a uniform 
solution to a traveling wave solution, one can readily rec- 
ognize the correspondence between our discrete simula- 
tions and corresponding two-fluid model calculations of 
Glasser et al. 0. 

Our discrete simulations clearly show that the solid 
phase stresses and the granular temperature obtained 
by ensemble-averaging over many realizations are path- 
dependent. We have checked that the normal stress dif- 
ference in bubbling fluidized beds is different from the 
well-known Burnett order effect in sheared granular ffuids 
(See Appendix). The traveling wave can be partitioned 
into two regions, one where the assembly undergoes dila- 
tion and the other where it undergoes compaction. The 
solid phase stresses and the granular temperature in these 
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two regions differ markedly (see Fig. Wfiile one does 
expect some difference in tlie beliavior observed in the 
two regions, our simulations reveal that the difference 
is much larger than what can be attributed to the bulk 
viscosity correction in the widely used kinetic theory for 
granular materials (KTGM). The kinetic theory is lim- 
ited to situations where the volumetric strain rate is small 
compared to the rate at which local microstructure equi- 
librates, so that the effect of the volumetric strain rate 
can be treated as a small perturbation from the equilib- 
rium microstructure. Our computations show that, in 
this simple traveling wave, the volumetric strain rate is 
not small (see Fig. I12|) . Although these traveling waves 
are not observed in real gas-fluidized beds, it is such one- 
dimensional traveling waves that give way to bubble-like 
voids, and so one can readily expect that the same effects 
will be observed in bubble-like voids as well. 

Large (scaled) volumetric strain rates can be expected 
to lead to nonlinear effects, and this is indeed seen. Both 
bulk and shear viscosities in the solid phase are found to 
depend on the volumetric strain rate (see Figs. ^] and 
I13|l , and such dependence is not captured by the KTGM 
(which only considers the case of small-scale volumet- 
ric strain rates). When two- fluid models fail to account 
for such nonlinear effects, they cannot be expected to 
reproduce the macroscopic properties of the void region 
accurately. It is well established that bubble-like voids 
readily form in two-fluid model calculations even with a 
simple phenomenological closure for the particle phase 
stress 0,0, 113; however, the shape of such voids is not 
consistent with those observed experimentally. A part 
of this problem stems from numerical errors associated 
with discretization of the two-fluid model equations 53| ; 
however, numerical accuracy is not the sole factor. In 
the dense emulsion phase in the immediate vicinity of 
the bubble-like void, one can expect large-scale volumet- 
ric strain rates (just as in the ID-TW examined in this 
study) and the failure to account for the accompanying 
nonlinear effects will impact the characteristics of the 
bubble-like voids. 

In many technological applications involving gas- 
particle flows, particles in the 50 - 100 /im size range 
are employed. Such particles, belonging to Geldart type 
A [s^, manifest modest levels of inter-particle cohesion 
and the relevance of this cohesive interaction on the flow 
characteristics has been debated extensively in the litera- 
ture (for example, see a review article by Sundaresan |^); 
however, concrete conclusions are yet to emerge. In the 
present study, we have also examined the effect of inter- 
particle cohesion on the structure of the traveling wave. 
The inter-particle cohesive force model itself is quite sim- 
ple, but it suffices to explore the generic effects. It is 
found that all of the effects observed in non-cohesive sys- 
tems persist when a modest level of cohesion is added 
(see Figs. [Island [T7|) . Cohesion lowers the particle vol- 
ume fraction in the plateau region and hence the am- 
plitude of the wave. It also slows down the speed of 
the wave propagation (see caption of Fig. I15|) - which 



(in multi-dimensional systems) affords greater exchange 
of gas between the void and the surrounding emulsion. 
Thus the beneficial effect of modest levels of cohesion can 
be understood. 

Although we have presented results only for Bo up to 
2, we found that wave solutions could be obtained even 
for larger Bo values. At larger Bo values (e.g.. Bo — 
4), the waves became irregular and when Bo >~ 8, the 
entire assembly traveled as a single plug. This is not a 
surprising result, as it is well known that highly cohesive 
particles (Geldart type C [s^) do not fluidize well. Thus 
the range of Bo over which beds can be fluidized without 
additional excitation (such as mechanical vibrations) is 
Bo <^ 8. Since the wave propagation is irregular for 
Bo >^ 3, the narrow range of < Bo <~ 3 is more 
appropriate for smooth fluidization and hence quantita- 
tive analysis through ensemble-averaging. At such levels 
of cohesion, the particle assembly is clearly in a state of 
tension when it undergoes dilation (see Figs. I15land ll7|l . 
It also clearly speaks for the need to include in two-fluid 
models a path-dependent model for the granular pres- 
sure. 
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Appendix: The normal stress differences in sheared 
gas-fluidized beds 

Normal stress difference can arise due to Burnett order 
effects (Sela and Goldhirsh 13 )• Theoretical prediction 
for normal stress difference attributable to Burnett or- 
der effects is available only for the dilute limit, and a 
validated Burnett order hydrodynamical theory for bub- 
bling fluidized beds (i.e. dense suspensions) is currently 
unavailable. We flrst estimated the normal stress differ- 
ences ((Tzz — axx) using a general relationship provided 
by Sela and Goldhirsh [33| (even though it was developed 
for the dilute limit), and then considered a couple of rel- 
evant computational experiments in order to probe if the 
stress differences that we observed in our fluidized beds 
were consistent with Burnett order effects. 

According to a general Burnett order expression for 
the stresses (Eq. (46) in Sela and Goldhirsh 34]), the 
normal stress difference in ID-TW in fluidized beds de- 
pends on (products of) gradients of solid phase contin- 
uum variables, such as duz/dz, dT/dz, and d'^T/dz'^. We 
determined that, for the waves in our simulations, the 
dominant contribution came from a term of the form 
C(j)X'^{duz/dz)'^, where C is a constant of order of unity 
and A is the mean free path. In dense suspensions, A is 
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FIG. 18: Solid phase continuum variables of beds of ideal 
particles (e = 1.0; n = 0.0; Bo = 0; St = 55; (jjavg = 0.44). 

much smaller than unity. We found this term to be much 
smaller than the (cr^z — axx) values that we obtained in 
our simulations (see, for instance, Fig.|Sl(b)). Thus, there 
is no basis for attributing the stress difference seen in our 
simulations to Burnett order effects. 

To test this issue further, we also considered two addi- 
tional computational experiments: 

(i) The Burnett order theory by Sela and Gold- 
hirsh "34*1 for homogeneously sheared granular flows 
in the dilute limit showed that the normal stress 
differences still persist in the elastic limit. This 
theory also predicts that the first normal stress 
difference ((ctzz — o'a;a;)/(Tr(£)/_D), in our case) is 

a rapidly increasing function of inelasticity. We 
first computed different diagonal components of the 
stress in usual fluidized beds, using ideal particles 
(e — 1.0, /i = 0.0) used in Fig. O and examined 
their dependence on the inelasticity. 

(ii) Using dissipative particles (e = 0.9, /i = 0.1), we 
applied homogeneous shear in the ?/— direction to a 
fluidized bed; if there were any systematic normal 
stress difference between the xx— and yy— compo- 
nents, it would be mostly attributed to the Burnett 
order effects (and/or collisional anisotropy ||60|). 

In a bed of ideal particles (case (i) above, shown in 
Fig. ll8|) . a noticeable normal stress difference (ctzz — o'xx) 
persists at a level comparable to the dissipative particle 
system. Comparing the results for the ideal and dissipa- 
tive systems, we saw that the ratio {ozz — '^xx)/ (Tr(£) / D) 

was not a strong function of inelasticity. This further 
confirms that Burnett order effect is unimportant. 

When we applied shear to a fluidized bed (case (ii), 
shown in Fig.ll9|l in the y— direction, we still saw signifi- 
cant differences between xx— (or yy—) and zz— compo- 
nent normal stresses, but those between xx— and yy— 
components were relatively very small. The maximum 
volumetric strain rate (duz/dz ~ 0.2) in Fig. 1191 is com- 
parable to the shear strain rate applied {dux/dy = 0.2). 



FIG. 19: Solid phase continuum variables in sheared fluidized 
beds (computed by employing Lees-Edwards boundary condi- 
tion in the y— direction). The shear strain rate dux/dy of 0.2 
is applied, (e = 0.9; y, = 0.1; Bo = 0; St = 55; (jjavg = 0.44). 

When we increased the shear strain rate to ^ 1.0, the 
bed became homogenized and the region of low particle 
volume fraction was not visibly recognizable any more. 
Sheared gas fluidized beds are known to fluidize homoge- 
neously when the shear rate reaches some critical value, 
depending on the superficial gas flow rate ^| . Our com- 
putation yielded comparable values for the critical shear 
rate, however, quantitative comparison was not possible, 
as the critical rate also depends on the bed height 6IJ. 

To summarize, we used a Burnett order expression p34| 
for normal stress differences (which was developed for the 
dilute limit) to estimate its prediction for dense suspen- 
sions (as this is the only available theory). We found 
that such prediction significantly underestimated normal 
stress difference seen in our wave simulations. Further- 
more, when we assessed the first normal stress difference 
in beds of perfectly elastic particles, we did not observe 
any noticeable decrease (from that in beds of inelastic 
particles). In order to computationally obtain the nor- 
mal stress difference in sheared dense suspensions (as op- 
posed to dilute granular flows in the absence of gas), 
we considered a sheared fluidized bed. We found that 
the normal stress difference between the xx— and yy— 
components (presumably attributed to Burnett order ef- 
fects and/or collisional anisotropy) is still much smaller 
than the difference between the xx— and zz— compo- 
nents which arises from compaction and dilation. Based 
on all these results, we conclude that {azz — <^xx) seen 
in our traveling wave simulations is not due to Burnett 
order effects or collisional anisotropy. 
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